####Overview####
# File name:    "replication_hope_et_al"
# Last Change:  28.11.2025
# Purpose:      Replication for "The Hidden Cost of Tax Regressivity at the Top" (BJPS)

####1. System setups#####
#Load Packages
# install packages from CRAN
p_needed <- c("sf","margins","readstata13", "egg","ggpubr","giscoR","interplot", "installr", "sf", "tidyverse", "stargazer", "texreg","varhandle","ggthemes","grid","cowplot","foreign","magrittr","haven","xtable","ggrepel","ggplot2",
              "Matrix","gridExtra","grid","broom","multcomp","interactions", "mediation","interflex", "AER", "ggridges", "ggcorrplot", "plotly","readxl","xlsx","scales","RColorBrewer")
packages <- rownames(installed.packages())
p_to_install <- p_needed[!(p_needed %in% packages)]
if (length(p_to_install) > 0) {
  install.packages(p_to_install)
}
lapply(p_needed, require, character.only = TRUE)

#Set Seed
set.seed(2789)

#Set working directory to current folder
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#Load Data
main_data <- readRDS(file = "replication_data.Rdata") 


####Table 1 -- Treatment Effect on Tax Policy Preference#####

summary(top_1 <- glm(P1_1 ~ factor(treat), data = main_data))  

#Middle
summary(middle_1 <- glm(P2_1 ~ factor(treat), data = main_data))  

#Bottom
summary(bottom_1 <- glm(P3_1 ~ factor(treat), data = main_data))  

#WTP
summary(wtp_1 <- glm(P4_1 ~ factor(treat), data = main_data))  

screenreg(l = list(top_1, middle_1, bottom_1, wtp_1), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Figure 2 -- Effect of Treatment on Perceptions of Tax Regressivity at the Top and Fairness of the Tax System####

#Plot Mechanism
summary(manip_1 <- glm(rich_higher_tax ~ factor(treat), data = main_data))

model_summary <- summary(manip_1)

coef_val <- model_summary$coefficients[2, 1]  # Coefficient
p_val <- model_summary$coefficients[2, 4]     # p-value

# Create significance stars
stars <- if (p_val < 0.001) {
  "***"
} else if (p_val < 0.01) {
  "**"
} else if (p_val < 0.05) {
  "*"
} else if (p_val < 0.1) {
  "."
} else {
  ""
}

# Format label
label <- paste0("Coef = ", round(coef_val, 2), " ", stars)

# Determine y-position for the bracket
y_max <- max(main_data$rich_higher_tax, na.rm = TRUE)
bracket_y <- y_max -1.2       # bracket
label_y <- bracket_y +0.3    # label above bracket

# Plot
manip_plot <- main_data %>%
  filter(!is.na(treat)) %>%
  ggplot(aes(x = treat, y = rich_higher_tax, color = treat)) +
  geom_jitter(width = 0.1, alpha = 0.1) +
  stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, color = "black") +
  annotate("segment", x = 1, xend = 2, y = bracket_y, yend = bracket_y) +  # horizontal line
  annotate("segment", x = 1, xend = 1, y = bracket_y, yend = bracket_y - 0.1) +  # vertical tick left
  annotate("segment", x = 2, xend = 2, y = bracket_y, yend = bracket_y - 0.1) +  # vertical tick right
  annotate("text", x = 1.5, y = label_y, label = label, size = 5) +
  scale_color_manual(values = c("#1b9e77", "#7570b3")) +
  theme_bw() +
  labs(
    x = " ",
    y = "Agree That Richer Americans\n Pay Higher Taxes"
  ) +
  theme(
    legend.position = "none",  # Or "top" if preferred
    text = element_text(size = 16),
    axis.title = element_text(face = "bold")
  )


#Plot Fairness Mechanism
summary(mech_1 <- glm(tax_fair ~ factor(treat), data = main_data))  

model_summary <- summary(mech_1)

coef_val <- model_summary$coefficients[2, 1]  # Coefficient
p_val <- model_summary$coefficients[2, 4]     # p-value

# Create significance stars
stars <- if (p_val < 0.001) {
  "***"
} else if (p_val < 0.01) {
  "**"
} else if (p_val < 0.05) {
  "*"
} else if (p_val < 0.1) {
  "."
} else {
  ""
}

# Format label
label <- paste0("Coef = ", round(coef_val, 2), " ", stars)

# Determine y-position for the bracket
y_max <- max(main_data$tax_fair, na.rm = TRUE)
bracket_y <- y_max -2.2       # bracket
label_y <- bracket_y +0.3    # label above bracket

# Plot
fair_tax <- main_data %>%
  filter(!is.na(treat)) %>%
  ggplot(aes(x = treat, y = tax_fair, color = treat)) +
  geom_jitter(width = 0.1, alpha = 0.1) +
  stat_summary(fun = mean, geom = "point", size = 3, color = "black") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width = 0.2, color = "black") +
  annotate("segment", x = 1, xend = 2, y = bracket_y, yend = bracket_y) +  # horizontal line
  annotate("segment", x = 1, xend = 1, y = bracket_y, yend = bracket_y - 0.1) +  # vertical tick left
  annotate("segment", x = 2, xend = 2, y = bracket_y, yend = bracket_y - 0.1) +  # vertical tick right
  annotate("text", x = 1.5, y = label_y, label = label, size = 5) +
  scale_color_manual(values = c("#1b9e77", "#7570b3")) +
  theme_bw() +
  labs(
    x = " ",
    y = "Agree US Tax System \n Is Fair"
  ) +
  theme(
    legend.position = "none",  # Or "top" if preferred
    text = element_text(size = 16),
    axis.title = element_text(face = "bold")
  )

#Combine
plot_grid(manip_plot, fair_tax, nrow = 1)

####Figure 3 -- Mediation Analysis of Treatment Effect on Preferences for Taxing the Middle Class####
#Drop NAs
main_data_full <- main_data %>%
  filter(!(is.na(tax_fair)),
         !(is.na(P2_1)))

# Step 1: Model the mediator
model.M <- lm(tax_fair ~ treat_bin, data = main_data_full)

# Step 2: Model the outcome
model.Y <- lm(P2_1 ~ treat_bin + tax_fair, data = main_data_full)

# Step 3: Mediation analysis
med.out <- mediate(model.M, model.Y, treat = "treat_bin", mediator = "tax_fair", boot = TRUE, sims = 1000)

# Summary
summary(med.out)

effect_df <- data.frame(
  Effect = c("ACME (Indirect)", "ADE (Direct)", "Total Effect"),
  Estimate = c(med.out$d0, med.out$z0, med.out$tau.coef),
  Lower = c(med.out$d0.ci[1], med.out$z0.ci[1], med.out$tau.ci[1]),
  Upper = c(med.out$d0.ci[2], med.out$z0.ci[2], med.out$tau.ci[2])
)

# Plot
ggplot(effect_df, aes(y = Effect, x = Estimate)) +
  geom_point(size = 3, color = "#7570b3") +
  geom_errorbarh(aes(xmin = Lower, xmax = Upper), height = 0, color = "#7570b3") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  theme_bw()  +
  labs( x = "Effect Size",
        y = NULL
  ) + theme(
    text = element_text(size = 18),
    axis.title = element_text(face = "bold")
  )

####Table C1 -- Balance Checks – Determinants of Receiving Treatment####
summary(balance <- glm(treat_bin ~ income + age + gender + college + employment + rightist + trust + party_affiliation 
                       , data = main_data))
screenreg(l = list(balance), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE)

####Figure C1 -- Distribution of Support for Tax Reforms for Highest, Middle, and Lowest Income Earners, By Treatment####

bottom_distr <- ggplot(main_data, aes(x = factor(P3_1), fill = factor(treat))) +
  geom_bar(position = "dodge") +
  labs(x = "Lower - Raise Taxes on Lowest Income Earners", y = "Count", fill = "Treatment Group") +
  theme_bw()
middle_distr <- ggplot(main_data, aes(x = factor(P2_1), fill = factor(treat))) +
  geom_bar(position = "dodge") +
  labs(x = "Lower - Raise Taxes on Middle Income Earners", y = "Count", fill = "Treatment Group") +
  theme_bw()
top_distr <- ggplot(main_data, aes(x = factor(P1_1), fill = factor(treat))) +
  geom_bar(position = "dodge") +
  labs(x = "Lower - Raise Taxes on Highest Income Earners", y = "Count", fill = "Treatment Group") +
  theme_bw()

#Combine
plot_grid(top_distr, middle_distr, bottom_distr, nrow = 3)

####Table C2 -- Treatment Effect on Satisfaction With Democracy####
summary(dem_1 <- glm(P5_1 ~ factor(treat), data = main_data))  

screenreg(l = list(dem_1), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))
####Table D1 -- Treatment Effect on Tax Policy Preferences Including Covariates####

summary(top_2 <- glm(P1_1 ~ factor(treat)+
                       income  + age + gender  + college + employment + rightist + trust + party_affiliation 
                     , data = main_data))
summary(middle_2 <- glm(P2_1 ~ factor(treat)+
                          income  + age + gender  + college + employment + rightist + trust + party_affiliation
                        , data = main_data))
summary(bottom_2 <- glm(P3_1 ~ factor(treat)+
                          income  + age + gender  + college + employment + rightist + trust + party_affiliation
                        , data = main_data))
summary(wtp_2 <- glm(P4_1 ~ factor(treat)+
                       income  + age + gender  + college + employment + rightist + trust + party_affiliation 
                     , data = main_data))  

screenreg(l = list(top_2, middle_2, bottom_2, wtp_2), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Table D2 -- Treatment Effect on Tax Policy Preferences, Speeders Excluded####
main_data_speeders <- main_data %>%
  filter(Duration >150)


summary(top_3 <- glm(P1_1 ~ factor(treat)
                     , data = main_data_speeders))
summary(middle_3 <- glm(P2_1 ~ factor(treat)
                        , data = main_data_speeders))
summary(bottom_3 <- glm(P3_1 ~ factor(treat)
                        , data = main_data_speeders))
summary(wtp_3 <- glm(P4_1 ~ factor(treat)
                     , data = main_data_speeders))  
screenreg(l = list(top_3, middle_3, bottom_3, wtp_3), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE,  
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Table D3 -- Treatment Effect Tax Policy Preferences, Laggards Excluded####
main_data_laggards <- main_data %>%
  filter(Duration <= 900)


summary(top_4 <- glm(P1_1 ~ factor(treat)
                     , data = main_data_laggards))
summary(middle_4 <- glm(P2_1 ~ factor(treat)
                        , data = main_data_laggards))
summary(bottom_4 <- glm(P3_1 ~ factor(treat)
                        , data = main_data_laggards))
summary(wtp_4 <- glm(P4_1 ~ factor(treat)
                     , data = main_data_laggards))  
screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Table D4 -- Treatment Effect on Answering ’Don’t Know'####
main_data_dk <- main_data %>%
  mutate(
    tax_fair_binary = case_when(
      is.na(tax_fair) ~ 1,
      T ~ 0
    ),
    prog_binary = case_when(
      is.na(prog) ~ 1,
      T ~ 0
    ),
    P1_1_binary = case_when(
      is.na(P1_1) ~ 1,
      T ~ 0
    ),
    P2_1_binary = case_when(
      is.na(P2_1) ~ 1,
      T ~ 0
    ),
    P3_1_binary = case_when(
      is.na(P3_1) ~ 1,
      T ~ 0
    ),
    P4_1_binary = case_when(
      is.na(P4_1) ~ 1,
      T ~ 0
    ),
    P5_1_binary = case_when(
      is.na(P5_1) ~ 1,
      T ~ 0
    )
  )
summary(top_5 <- glm(P1_1_binary ~ factor(treat)
                     , data = main_data_dk))
summary(middle_5 <- glm(P2_1_binary ~ factor(treat)
                        , data = main_data_dk))
summary(bottom_5 <- glm(P3_1_binary ~ factor(treat)
                        , data = main_data_dk))
summary(wtp_5 <- glm(P4_1_binary ~ factor(treat)
                     , data = main_data_dk))  

screenreg(l = list(top_5, middle_5, bottom_5, wtp_5), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))


####Table D5 -- Treatment Effect on Tax Policy Preferences, Adjustment for Multiple Comparisons####
#Top
summary(top_1 <- glm(P1_1 ~ factor(treat), data = main_data))  

#Middle
summary(middle_1 <- glm(P2_1 ~ factor(treat), data = main_data))  

#Bottom
summary(bottom_1 <- glm(P3_1 ~ factor(treat), data = main_data))  

#WTP
summary(wtp_1 <- glm(P4_1 ~ factor(treat), data = main_data))  

pvals <- c(
  coef(summary(top_1))[, "Pr(>|t|)"][-1],   
  coef(summary(middle_1))[, "Pr(>|t|)"][-1],
  coef(summary(bottom_1))[, "Pr(>|t|)"][-1],
  coef(summary(wtp_1))[, "Pr(>|t|)"][-1]
)

# Bonferroni correction
pvals_bonf <- p.adjust(pvals, method = "BH")
top_texreg <- createTexreg(
  coef.names = rownames(coef(summary(top_1))),
  coef = coef(top_1),
  se = coef(summary(top_1))[, "Std. Error"],
  pvalues = c(coef(summary(top_1))[1, "Pr(>|t|)"], pvals_bonf[1])
)

middle_texreg <- createTexreg(
  coef.names = rownames(coef(summary(middle_1))),
  coef = coef(middle_1),
  se = coef(summary(middle_1))[, "Std. Error"],
  pvalues = c(coef(summary(middle_1))[1, "Pr(>|t|)"], pvals_bonf[2])
)

bottom_texreg <- createTexreg(
  coef.names = rownames(coef(summary(bottom_1))),
  coef = coef(bottom_1),
  se = coef(summary(bottom_1))[, "Std. Error"],
  pvalues = c(coef(summary(bottom_1))[1, "Pr(>|t|)"], pvals_bonf[3])
)

wtp_texreg <- createTexreg(
  coef.names = rownames(coef(summary(wtp_1))),
  coef = coef(wtp_1),
  se = coef(summary(wtp_1))[, "Std. Error"],
  pvalues = c(coef(summary(wtp_1))[1, "Pr(>|t|)"], pvals_bonf[4])
)



screenreg(l = list(top_texreg, middle_texreg, bottom_texreg, wtp_texreg), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Table D6 -- Treatment Effect on Support for Raising Taxes on Lowest Income Earners, Dichotomised####
main_data_bottom <-
  main_data %>%
  mutate(bottom_bin = case_when(
    P3_1 == 0 ~ 0,
    P3_1 > 0  ~ 1,
    T ~ P3_1
  ))

summary(bottom_1 <- glm(bottom_bin ~ factor(treat), data = main_data_bottom))  

screenreg(l = list(bottom_1), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment"))

####Table E1 -- Treatment Effect on Tax Policy Preferences, Interaction Gender####
summary(top_4 <- glm(P1_1 ~ factor(treat) * gender
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * gender
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * gender
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * gender
                     , data = main_data))  
screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "genderMale"  = "Male (Ref. Female)", 
                              "genderOther"  = "Other (Ref. Female)" , 
                              "genderPrefer not to answer"  = "Prefer not to answer (Ref. Female)",
                              "factor(treat)Treatment:genderMale"  = "Treatment*Male (Ref. Female)", 
                              "factor(treat)Treatment:genderOther"  = "Treatment*Other (Ref. Female)" , 
                              "factor(treat)Treatment:genderPrefer not to answer"  = "Treatment*Prefer not to answer (Ref. Female)"))

####Table E2 -- Treatment Effect on Tax Policy Preferences, Interaction Age####

summary(top_4 <- glm(P1_1 ~ factor(treat) * age
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * age
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * age
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * age
                     , data = main_data))  
screenreg(l = list(top_4, middle_4, bottom_4,wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "age"  = "Age", 
                              "factor(treat)Treatment:age"  = "Treatment*Age"))

####Table E3 -- Treatment Effect on Tax Policy Preferences, Interaction Income####

summary(top_4 <- glm(P1_1 ~ factor(treat) * income
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * income
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * income
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * income
                     , data = main_data))
screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "income"  = "Income", 
                              "factor(treat)Treatment:income"  = "Treatment*Income"))

####Figure E1 -- Treatment Effect on Tax Policy Preferences by Income, Binning Estimator####
main_data <- data.frame(main_data)

rich_inter <- plot.interflex(interflex(Y = "P1_1", D = "treat_bin", X = "income", 
                                       data = main_data, estimator = "binning",
                                       na.rm = T, nbins = 3),
                             ylab = "Marginal Effect of Treatment", 
                             xlab ="Income Groups", theme.bw = T,
                             main = "Tax Top",
                             text = element_text(size = 12))


middle_inter <- plot.interflex(interflex(Y = "P2_1", D = "treat_bin", X = "income", 
                                         data = main_data, estimator = "binning",
                                         na.rm = T, nbins = 3),
                               ylab = "Marginal Effect of Treatment", 
                               xlab ="Income Groups", theme.bw = T,
                               main = "Tax Middle",
                               text = element_text(size = 12))

bottom_inter <- plot.interflex(interflex(Y = "P3_1", D = "treat_bin", X = "income", 
                                         data = main_data, estimator = "binning",
                                         na.rm = T, nbins = 3),
                               ylab = "Marginal Effect of Treatment", 
                               xlab ="Income Groups", theme.bw = T,
                               main = "Tax Bottom",
                               text = element_text(size = 12))

wtp_inter <- plot.interflex(interflex(Y = "P4_1", D = "treat_bin", X = "income", 
                                      data = main_data, estimator = "binning",
                                      na.rm = T, nbins = 3),
                            ylab = "Marginal Effect of Treatment", 
                            xlab ="Income Groups", theme.bw = T,
                            main = "WTP",
                            text = element_text(size = 12))

plot_grid(rich_inter, middle_inter, bottom_inter, wtp_inter, nrow = 1)
####Table E4 -- Treatment Effect on Tax Policy Preferences, Interaction College####

summary(top_4 <- glm(P1_1 ~ factor(treat) * college
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * college
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * college
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * college
                     , data = main_data))  

screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "college"  = "College", 
                              "factor(treat)Treatment:college"  = "Treatment*College"))
####Table E5 -- Treatment Effect on Tax Policy Preferences, Interaction Left-Right####


summary(top_4 <- glm(P1_1 ~ factor(treat) * rightist
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * rightist
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * rightist
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * rightist
                     , data = main_data))  

screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "rightist"  = "Left-Right", 
                              "factor(treat)Treatment:rightist"  = "Treatment*Left-Right"))
####Table E6 -- Treatment Effect on Tax Policy Preferences, Interaction Party Affiliation####
summary(top_4 <- glm(P1_1 ~ factor(treat) * party_affiliation
                     , data = main_data))
summary(middle_4 <- glm(P2_1 ~ factor(treat) * party_affiliation
                        , data = main_data))
summary(bottom_4 <- glm(P3_1 ~ factor(treat) * party_affiliation
                        , data = main_data))
summary(wtp_4 <- glm(P4_1 ~ factor(treat) * party_affiliation
                     , data = main_data))  
screenreg(l = list(top_4, middle_4, bottom_4, wtp_4), include.ci = FALSE,
       stars = c(0.01, 0.05, 0.1), include.bic = TRUE, include.aic = TRUE, include.rsquared = F, include.variance = FALSE,  include.deviance = F,
       booktabs = TRUE, dcolumn = TRUE, longtable = TRUE, center = TRUE, digits = 4, include.missings = FALSE, no.margin = TRUE, 
       custom.coef.map = list("factor(treat)Treatment" = "Treatment",
                              "party_affiliationRepublican Party"  = "Republican (Ref. Democrat)", 
                              "party_affiliationOther"  = "Other (Ref. Democrat)", 
                              "factor(treat)Treatment:party_affiliationRepublican Party"  = "Treat.*Republican (Re. Democrat)",
                              "factor(treat)Treatment:party_affiliationOther"  = "Treat.*Other (Ref. Democrat)"))












